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An efficient Monte Carlo method is extended to evaluate directly domain-wall free-energy for ran- 
domly frustrated spin systems. Using the method, critical phenomena of spin-glass phase transition 
is investigated in 4d± J Ising model under the replica boundary condition. Our values of the critical 
temperature and exponent, obtained by finite-size scaling, are in good agreement with those of the 
standard MC and the series expansion studies. In addition, two exponents, the stiffness exponent 
and the fractal dimension of the domain wall, which characterize the ordered phase, are obtained. 
The latter value is larger than d — 1, indicating that the domain wall is really rough in the 4d Ising 
spin glass phase. 



PACS numbers: 75.50.Lk,75.40.Mg 
I. INTRODUCTION 

Numerical simulations, in particular Monte Carlo 
(MC) methods, have played a quite important role in 
spin glass (SG) studies! For example, very large-scale 
MC simulations have strongly suggested the existence 
of a SGfpliase transition in three-dimensional Ising SG 
systems an In these studies, a cumulant of SG over- 
lap function q, so called the Binder parameter, has fre- 
quently used in order to extract critical temperature T c . 



However, the Bin 
(EA) Ising model 



a parameter in 3d Edwards- Anderson 
merely depends cm system sizes be- 
low T c as compared to that above T C Q. Moreover, un- 
usual size dependence of the Binder parameter is ob- 
served in ^short-ranged Ising EA model under the mag- 
netic fieldul. Consequently, the existence of the SG phase 
transition under the field has still remained unclear. In 
order to settle the issue and make a progress toward well- 
understanding of SG picture, we consider other numerical 
analyzes to be quite necessary. 

In this direction, some promising ways have recently 
been proposed based on non-equilibrium dynamicsaO and 
the idea of non-self-averagingl3, but here we pay at- 
tention to domain- wall renormalization group, method 
(DWRG) originally proposed by McMillaiM The 
DWRG estimates a singular part of free energy by calcu- 
lating the domain-wall free energy, AF, which is defined 
as the free-energy difference between periodic- and anti- 
periodic boundary conditions (BC). In the scaling regime 
at low temperatures, AF follows a power law as a func- 
tion of the system size L, AF ~ L e , where the stiffness 
exponent 9 is related to the rigidity of the system. If the 
exponent 9 takes a positive value at a temperature, then 
the system stays in an ordered phase. On the other hand, 
a negative exponent means a disordered phase. In this 
sense, the sign of the exponent 9 is an indicator of the 
existence of long range ordering. This exponent 9 also 
characterizes a low-energy excitation in the SG phase 
and is predicted to be smaller than (d — llZ^c? being 
dimensionality, in the droplet scaling theory 12 



The DWRG approach relies on an accurate way for 
estimating the free-energy difference between two BCs. 
It is a difficult task in general for MC method to esti- 
mate free energy or entropy. Except for the numerical 
transfer matrix method for Ising models, it is therefore 
usual to integrate over the free-energy derivative, mea- 
sured by MC simulations, along a parameter path be- 
tween a reference system and the one of interest. As for 
zero-temperature calculations, various optimization tecit, 
niques have been demonstrated to be useful for IsingJiJ'E£l 
and vector spin systemsMI. These facts restrict so far to 
rather small sizes and/or at zero temperature. In this 
work we have developed— a boundary-flip MC method 
proposed by HasenbuschEZI which allows us to estimate 
free-energy difference at a finite temperature directly from 
MC simulation. In applying a naive boundary-flip MC 
method to large systems and/or at low temperatures, one 
may encounter a hardly relaxing problem even in sim- 
ple models without many meta-stable states, namely the 
system is trapped-into a local area in the phase space. 
The original worktU has successfully overcome the relax- 
ational problem by combining the method with the clus- 
ter MC dynamics. 

In the present paper, we have proposed an alterna- 
tive strategy, which is the boundary-flip method with 
exchange MC (EMC) methodla, in order to make the 
relaxation faster. This combined method is found to be 
quit efficient for randomly frustrated spin systems such as 
spin glasses, while the original method based on the clus- 
ter MC method is restricted to non- frustrated systems. 
The present method is applicable to a wide class of spin 
systems. Moreover, the direct measurements have an 
advantage over the thermodynamic integration method 
from a numerical standpoint, because statistical error 
is controlled within MC scheme in the former. Conse- 
quently, we have succeeded to estimate the free-energy 
difference in a SG model, accurately enough to observe 
systematic correction to finite-size scaling. 

For applying DWRG to SG systems, we need to 
choose the relevant boundary conditions to the ordered 
phase. The standard approach has often used a ran- 
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domly fixed spin boundary conditionE3. Instead, we ena.- 
ploy the replica boundary condition proposed by OzekicJ, 
in which two real replicas are coupled with each other 
through a boundary surface. The replica boundary con- 
dition provides that the domain-wall free energy becomes 
positive at any bond disorder, implying that it is conju- 
gates to the SG ordering. This positivity is of benefit to 
us for estimating the domain-wall free energy accurately 
from a numerical point of view. 

Here we study the Ad ± J Ising SG model under the 
replica BC by the novel MC method. We obtain the crit- 
ical temperature and the exponent by finite-size-scaling 
analysis of the domain-wall free energy, in agreement 
with the previous works. In addition, we estimate two 
exponents, the stiffness exponent 6 and the fractal di- 
mension d s of the domain wall. We find that d s is larger 
than d — 1 in the SG phase. 

This paper is organized as follows: In the next section, 
we explain the method for calculating the domain-wall 
free energy. Section |l| is mainly devoted to discussion 
about the replica boundary conditions. We give an inter- 
pretation of domain wall appearing in the replica bound- 
ary condition and propose a way to measure the morphol- 
ogy of the domain wall. We show results for application 
of the method to 4d± J Ising SG model in the section IV. 
In the last section, possible extensions of the method and 
nature of the low-temperature phase are discussed. Ap- 
pendix A contains a way for setting temperature points 
which is needed before simulation in the exchange MC 
method. 



II. BOUNDARY FLIP MC METHOD WITH 
EXCHANGE PROCESS 

In this section we describe a method that allows us 
to evaluate directly the domain-wall free-energy using 
MC simulations. For the sake of simplicity, we restrict 
ourselves to Ising spin systems and fixed spin boundary 
conditions. Let us consider a total model Hamiltonian 
defined by 

Wtot(f) Si, S 2 ) = W m odel(f) + «^bc(^ Si, S 2 ), (1) 

where a denotes Ising spin variable defined on a d- 
dimensional hyper-cubic lattice V of L d and two addi- 
tional spins Si and S 2 represent boundary spins. The 
second term gives a coupling between the model system 
and the boundary spins along one direction as 



Hbc(o~, Si,S 2 ) 



E 



Ji.l°~i Si 



E 

i£d 2 V 



Ji,2PiS 2 , (2) 



where the summation runs over one surface d{V of the 
lattice V and its opposite surface d 2 V. A standard peri- 
odic boundary condition for <j is used along the remaining 
directions. Then the total partition function Z to t and the 
free energy F to t of this whole system are defined by 



Ztot(T) = Tr^ Su s 2 }exp(-H tot (a, Si,S 2 )/T) 

= exp(-F tot (T)/T), (3) 

where T is temperature and we set the Boltzmann con- 
stant to unity. The phase space of the total Hamiltonian 
is enlarged by adding the degree of freedom of the bound- 
ary spins Si and S 2 . When these spins are in parallel, 
the boundary condition is regarded as periodic and sim- 
ilarly the anti-periodic boundary condition corresponds 
to anti-parallel boundary spins. For a given temperature 
the probability for finding the periodic boundary condi- 
tion is given by 



Pp(T) 



_ ^{a : S 1 ,S 2 } S S 1 ,S 2 exp(-W t ot(CT,Si, S 2 )/T) 



Z to t{T) 



Zp(T) 



(4) 



where S is the Kroneker delta function. This quantity is 
accessible from MC simulation, namely it is nothing but 
the probability for realizing the periodic BC during MC 
simulation in which the boundary spins as well as the 
bulk spins are updated according to a standard MC pro- 
cedure. In terms of the probability and the corresponding 
one to the antiperiodic BC, domain-wall free-energy AF 
we want to investigate is given by 



exp(/3AF(T)) 



3 -/3(F P -F AP ) 



Z P Pp(T) 
Z A p Pap CO' 



(5) 



This is the basic idea of-the boundary flip MC method 
proposed by HasenbuschliZI. When we adopt a naive local 
updating process for the boundary spins in the boundary- 
flip MC method, however, we are at once faced to a 
hardly relaxing problem. For example, once the anti- 
periodic boundary conditions and the domain-wall struc- 
ture in the system are realized in the simulation at low 
temperatures, as shown in Fig. ||, the boundary spins 
are kept to be fixed in the sense that the probability for 
flipping these spins is vanishing in practice. This fact 
makes statistical error of AF significantly large. The 
original workO has overcome this difficulty by utilizing 
the modified cluster flip. We can also practically solve 
this so-called hardly relaxing problem using recently pro- 
posed extended ensemble methods such as the jmilti- 
canonical MC methodo, the simulated tempering^ and 
the exchange MC methodlia. In fact, a similar difficulty 
has been overcome using .the multicanonical idea in the 
lattice-switch MC methocc-3, which has been proposed to 
estimate the free-energy difference between two different 
crystalline structure in a hard sphere system. 

In the present work, we employ the exchange MC 
method (EMC) in order to obtain an efficient path 
between two boundary-condition states. In the EMC 
method, we simulate a combined system which consists of 
non-interacting M replicated system. The m-th replica 
is simulated independently with its own external variable 
such as temperature. We introduce exchange process be- 
tween configurations of two of M replicas with the whole 
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FIG. 1. Typical example for metastable configuration in 
a ferromagnetic model. 



a being unity, as shown in Fig. |J. It is noted that the 
parameter region of the our final interest lies on the line 
with a = 1 around T c and below. An efficient choice 
of the exchange line would depend on systems we want 
to investigate. Actual implementation to the Ising spin 
glass model will be explained in detail in Section IV. 



III. REPLICA BOUNDARY CONDITIONS FOR 
SG SYSTEMS 
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FIG. 2. Schematic picture of the exchange line in param- 
eter space. 



combined system remaining in equilibrium. One possible 
way for obtaining the path is that we distribute various 
values of the coupling a in eq. ([!]) ranging from to 1 to 
M replicas. A target system we are physically interested 
in is the replica with a = 1. For a replica with null cou- 
pling of a, which we call a source system, the boundary 
spins can be flipped freely. Therefore, the path between 
different boundary condition states in the target system 
would be recovered by the exchange process through the 
source system. 

In randomly frustrated spin systems such as SG mod- 
els, there is another serious relaxation problem arising 
from bulk spins in the model system itself. Thds problem 
can be overcome also by the EMC method.lia When we 
distribute M temperature points widely including high 
enough temperature in a disordered phase, configurations 
at low temperatures are expected to be refreshed through 
the exchange process. The EMC method, has turned out 
to work efficiently in the SG systemaiflO Therefore, for 
the boundary- flip MC method on SG models, we need to 
construct the EMC method in two-dimensional parame- 
ter space of the coupling a and the temperature T. It 
is possible to introduce the exchange process in the two 
parameter space, but it is quite time consuming. In the 
present work, therefore, we choose an exchange line in 
the two dimensional space appropriately, namely we set 
a system at high temperature with a = as one end of 
the exchange line and systems at lower temperature with 



In this section, we discuss how to choose boundary 
condition for SG systems in the DWRG study. We con- 
centrate ourselves on a way of choice of boundary con- 
dition along one direction while the remaining ones are 
considered to be. given appropriately. In a conventional 
DWRG studiesEJa as well as the defect energy method, 
a boundary condition frequently used is a connected spin 
BC in which the corresponding boundary term in eq. (|l|) 
is described by 



Hbc = 

ied 1 v,jed 2 v 



(6) 



The case with a/\a\ = 1(— 1) is regarded as the (anti-) 
periodic boundary condition. For the boundary condi- 
tion defined by eq. (^|), the boundary-flip MC method 
can be applied by treating the sign of the coupling a as a 
MC dynamical variable. In SG systems, the free-energy 
difference between these BCs cannot be assured positive 
so that the width of distribution of the free-energy dif- 
ference is examined as an effective coupling of the SG 
ordering, F e g — y/ (Fap — Fp) 2 - To evaluate the mean 
width is rather difficult as compared with the average in 
numerical calculations. Further, it is less clear how the 
domain wall is created in a random spin system under 
these BCs. 

In order to avoid the difficulty and make clear an 
idea of the domain wall, Ozekin3 has proposed a replica 
boundary condition (RBC), in which two real replicas are 
prepared with the same bond realization. Its essential 
point is to introduce a uniform coupling between these 
two replicas only for one surface doV along a given direc- 
tion. For the other directions periodic BC is employed as 
usual. We show explicitly an example expressed as the 
Ising Hamiltonian 



udel 



for) 



^2 Jij(o-i<Tj 



E 



(ij) 



J int y ' T i 1 i ■ 

iGd V 



(7) 



where both a and r are Ising variables and the summa- 
tion of the first term runs over nearest neighbor bonds. 
The second term corresponds to the replica interaction 
mentioned above. When Jj nt is set to (anti-) ferromag- 
netic, the boundary condition is called replica (anti-) pe- 
riodic, R(A)PBC. Spins on the opposite side of 8qV are 
kept randomly fixed with Oi = Tj. 
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Ferromagnetic interactions between the replicas in the 
RPBC prefer a self-overlap state, even if the system has 
many local minima or pure states. Namely, one replica 
gives an effective conjugated field to the other replica 
through the inter-replica interaction, ft is convenient to 
consider the domain wall in terms of the replica overlap 
qi = <JiTi. The self-overlap state is characterized by posi- 
tive values of qi at all the sites, meaning no domain wall 
in the system. At sites on the opposite surfaces of 8qV, 
qi take unity by definition, irrespective of R(A)PBCs. 
On the other hand, anti-ferromagnetic inter-couplings be- 
tween the replicas in the RAPBC would induce negative 
overlap at sites near the coupling. Therefore, at least 
one domain wall, characterized by a region where the 
sign of qi changes, likely appears in the RAPBC, if the 
system has a rigid ordered state. From a mathematical 
point of view, non-negativity of the free-energy difference 
AFr = Frapbc — -Frpbc under the replica BC has been 
proven rigorously in any random Ising model at any-fi,- 
nite temperature using the transfer matrix formalisircJ. 
This non-negativity holds true irrespectively of a choice 
of spins on the surface opposite to d V. As a result, only 
average of the domain- wall free energy is needed for esti- 
mating an relevant effective coupling of the SG ordering. 
This is advantageous for reducing statistical error of AFr 
from which a transition point from paramagnetic to SG 
phase is detected. 

An additional merit of the replica boundary condition 
is that we can discuss the morphology of the domain wall 
at finite temperatures. In terms of the local overlap g,-, 
the area of the domain boundary mentioned above is ex- 



pressed as W = 



(ij) 2< 



where the summation 



is over nearest-neighboring pairs. . Then we can extract 
directly domain-wall properties such as its fractal dimen- 
sion, from the difference AW(T) defined by 

AW(T) = - y^((o-igjTjTj)RPBC - (0-i0-jTiTj) R APBc), 



(ij) 



(8) 



where (• • -)r(a)pbc denotes the thermal average under 
the replica (anti-) periodic BC. This quantity is also re- 
garded as a difference of link correlationt3 between two 
boundary conditions in ± J models. The correlation func- 
tion as well as the replicaj-averlap have been studied in a 
similar replicated systemE£l which has a global coupling 
between the replicas. This coupled system is different 
from the present system under RBC. In particular, the 
correlation function (^) is related to domain-wall prop- 
erties only in the RBC. The domain-wall area AW has 
not been directly studied so far in SG systems, except 
for the zero tenaBerature calculation in two dimensional 
Ising SG modelEi We will present new results for AW 
in the next section. 



IV. RESULTS 

In this section, we present results of an application of 
the novel MC method explained in the previous sections 
to Ad ± J Ising SG model. The interactions {Jij} in 
eq. (Q) are random variables which take values ± J with 
equal probability. The boundary-flip MC method can 
be applied to the replica BC by regarding the sign of 
the interaction Jj nt in eq. (0) as a dynamical variable. 
Equivalently these boundary conditions are defined by 
relative direction of the boundary spins Si and 5*2 added 
to eq. (0) whose J; n t are fixed to be positive. Then, the 
boundary part in eq. (El) is given by 



7~Ibc(°~, t, Si,S 2 ) 



iedV 



Ji{uiS\ + nS2), 



(9) 



where the interactions Jj are also distributed randomly. 
In the present work we adopt this method with the 
boundary spins. 

The number of the replicas M in the EMC method is 
fixed 32 irrespectively of the system sizes to utilize the 
multi-spin coding technique. Each replica with the pa- 
rameters a and T tries to exchange configuration with 
the nearest replica in the parameter space. As we have 
explained in Section ffS we choose in this two-parameter 
space, a line on which M replicas are prepared. The line 
chosen is such that value of a is unity below a certain 
temperature T m , but it decreases like a Gaussian formula 
as a function of T — T m above T m . The onset T m is set 
to be about 2 times of the critical temperature. We dis- 
tribute the set of the parameters to the 32 replicas such 
that the acceptance ratio for each exchange process be- 
comes independent of the replicas. This can be succeeded 
by a simple iteration method using the energy function, 
which is estimated from a short preliminary run. Details 
of the iteration method is explained in Appendix. 

As an equilibration check, we study time evolution 
of AFr starting from two initial conditions: periodic 
and anti-periodic boundary conditions imposed for the 
whole replicated systems in the EMC simulation. The 
initial conditions for the bulk spins are chosen at ran- 
dom. The free-energy difference AFr is estimated as a 
function of MC step t by averaging over short MC steps 
around the time t. In the case of the whole anti-periodic 
BC, free-energy difference, starting from a large negative 
value at the initial time, evolves toward equilibrium. The 
other estimation with the periodic BC at the initial time 
reaches to the equilibrium value from the opposite direc- 
tion to that of the anti-periodic BC. In equilibrium, two 
curves coincide with each other. As expected, we see in 



Fig. IV that the equilibration of AFr is obtained after 
a certain time. It should be noted that the relaxational 
function approaching to the equilibrium value follows an 
exponential law rather than a power law observed in the 
standard SG simulations. This implies the existence of a 
typical time scale for equilibration in the present method. 
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FIG. 3. The domain-wall free-energy of 4d ± J Ising SG 
model with L — 8 and T — 1.694 well below the SG transi- 
tion temperature as a function of MC steps. The upper data 
marked by open triangle are started from the periodic bound- 
ary condition for the whole system, while the lower one from 
the anti-periodic one. Each point at the time t is obtained 
by averaging over 2,000 MCS around t and error bars are 
estimated from statistical fluctuation over 10 samples. 



We thus expect that the system really reaches equilib- 
rium after a few times of such time scale. We estimated 
the time scale for other sizes and determined the MC 
steps (MCS) for thermalization and measurements. For 
example, in simulations of Ad case with L = 8, we take 
9.6 x 10 4 MCS for the initial step and 2.0 x 10 5 MCS 
for measurement. We have also checked that the ergodic 
timeEM is about 3 x 10 2 , 3.0 x 10 4 , 5.8 x 10 4 and 1.7 x 10 5 
MCS on average for L — 4,6, 8, and 10, respectively. 

We show temperature dependence of AFr for the 4d 
Ising SG model in Fig. IV. The lattice size studied are 
L = 4, 6, 8, and 10 with samples 2197, 2060, 1332, and 
892, respectively. According to the standard finite-size- 
scaling argument, the domain- wall free energy should be 
scaled as 



AF R (L,T)~F ((T-T c )L 1 /»), 



(10) 



where the parameter v denotes the critical exponent of 
the correlation length and Fo is a scaling function. There- 
fore, the critical temperature can be located in the point 
where AFr for different sizes as a function of T cross 
with each others. The crossing feature of AFr at T c is 
common to the Binder parameter. In fact, as shown in 
Fig. IV, crossing of AFr of two different sizes is seen at a 
certain temperature. However the crossing temperature 
is found to shift systematically to low temperature side 
as the system size increases, implying that correction to 
the finite-size scaling is significant. Wc consider correc- 
tion due to the leading irrelevant scaling variable whose 
scaling dimension is to; 

AF R (L, T) ~ F ((T - T C )L X I V ) + L— fl((T - T c )L l ' v ), 

(11) 
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FIG. 4. Temperature dependence of the domain-wall 
free-energy for 4d ±J Ising SG model near the critical tem- 
perature. These lines are for guide of eyes. 

These exponent v and u> and the critical temperature 
T c are determined by fitting the simulated data to the 
scaling formula (]Tl|), where the scaling functions Fq and 
Fi are assumed to be given by third order polynomial 
functions. From the fitting, we estimate T c — 2.00(4), 
v = 0.92(6), and u> — 1.5(9). The finite size scaling of 
Fo after subtraction of the leading correction is plotted 
in Fig IV, where all the data points are found to col- 
lapse almost into a universal function. The scaling plot 
including the smallest size L = 4 is obtained only when 
the leading term of the correction is taken into account. 
The estimated critical temperature is consistent, with the 
previous results obtained bp-JJae MC methodtll and the 
high temperature expansiorjcJ^. Our result for v is also 
in agreement with these expansion studies, and not vere 
different with that obtained, by MC simulations for ± Jell 
and Gaussian distribution^. Since the system sizes used 
in the present work are larger than those in the previous 
MC simulations, we expect that our estimation is reli- 
able. The irrelevant exponent u> is, to our knowledge, 
the first estimation for 4d Ising SG model by MC simu- 
lation, but its value is slightly lower than that obtained 
from the series expansiontJ that quoted about 3. 

At low enough temperature, the domain-wall free en- 
ergy is expected to be scaled as 



AFr(L, T) ~ L 6 



(12) 



where 9 is an exponent which gives the characteristic 
energy scale L e of low energy excitations of typical size 
L. We cannot evaluate AFr at low temperatures enough 
to distinguish the low temperature properties from the 
critical behavior. Here we try to estimate the exponent 
9 from the scaling function of AFr . We assume that the 
behavior of AFr at a large length scale is also described 
by the scaling form of cq. ( pr[ ) near below T c . This 
assumption implies that the asymptotic behavior of the 
scaling function F is predicted as 



F (x) 



(13) 
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FIG. 5. Finite-size scaling plot of the domain-wall free 
energy in 4d ± J Ising SG model. The leading correction 
to the scaling is taken into account. The scaling plot after 
subtraction of the leading correction is shown. The estimated 
scaling parameters are T c — 2.00(4), v = 0.92(6) and the irrel- 
evant exponent uj — 1.5(9). The slope of the scaling function 
is asymptotically close to 0.75(1), meaning that the stiffness 
exponent 9 is 0.82(6). 



at x — » — oo. We examine this scaling idea in the simple 
3c? Ising ferromagnetic model, where the stiffness expo- 
nent coincides with the surface dimensions d — 1. We 
estimate the domain-wall free energy by the present MC 
method under the connected spin BC described in (^). 
In the 3d Ising model, we scale the data to the lead- 
ing scaling formula ( jlO| ) without the correction, because 
we have not observed a shift of the crossing temperature 
under our numerical accuracy. The finite-size scaling of 
the domain-wall free energy works well as observed in 
Fig. [| The asymptotic behavior of the scaling function 
gives 9v ~ 1.27, compatible with the well-known values 
of v and 9 = d — 1. 

Let us turn to the 4c? Ising SG model. The stiffness 
exponent 9 in SG systems is expected much smaller than 
that of the ferromagnetic model. The droplet theory pre- 
dicted the upper bound of 9 to be (d— l)/2cj. We extract 
value of 9 from the scaling function obtained in Fig. 
We fit the scaled data with the scaling variable x larger 
than 3 to a power law. The best fit is obtained with the 
exponent 9v = 0.75(1), which yields the stiffness expo- 
nent of 9 = 0.82(6). 

We also investigate the domain-wall area AW defined 
by eq. (||) in this model, which is easily calculated in 
the present MC scheme. A scaling analysis similar to the 
one for AFr is performed for AW, taking into account 
the leading correction to the scaling. It is noted that in 
contrast with the AFr scaling, AW is proportional to 
L 2 I V near T c because it has essentially the same scaling 
dimension as the energy-energy correlation function. The 
finite-size-scaling plot for AW is shown in Fig. IV , where 
the critical temperature is used which is estimated by the 
AFr scaling. The scaling nicely works both above and 
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FIG. 6. Finite size scaling plot of the domain- wall free en- 
ergy in 3d ferromagnetic Ising model. The parameters of the 
scaling are estimated as follows: T c = 4.5117(4), v = 0.624(7). 
The asymptotic behavior of the scaling function follows a 
power law as a function of the scaling parameter (T — T c )L 1//y 
with slope 1.27(1). The value of the slope is compatible with 
the low temperature behavior, namely 9v, being 9 = d — 1. 



below T c and the estimated v value is consistent with 
that from AFr. We suppose that at low temperature 
the domain wall in the SG system is rather rough. Cor- 
respondingly the domain-wall area AW is expected to 
follow a power law on size with a non trivial fractal di- 
mension d s . We estimate d s by extracting the asymptotic 
behavior of the scaling function of AW in the same way 
as in the analysis of AFr. The asymptotic slope of the 
scaling function is d s v — 2. The fractal dimension of this 
model is found tCf-be 3.13(2). According to the Bray- 
Moore scaling lawtS, the exponents 9 and d s are related 
to the chaos exponent £ 



S 2 



(14) 



By this combined with the values of 9 and d s obtained 
here, our estimation of ( is 0.75(6). This value is sma 
than those of MC simulations for 4c? Ising SG modelstJ 
but rather close to that b« the Migdal-Kadanoff renor- 
malization group analysisEH. 



V. DISCUSSION AND SUMMARY 

We have developed a numerical method which enable 
us to estimate free-energy difference directly from MC 
simulation. It is a boundary flip MC method, in which 
the replica boundary conditions and the exchange MC 
technique are incorporated. The proposed method works 
well in the short-range Ising SG model. This method pre- 
sented here can be applied to various spin systems includ- 
ing vector spin models because our argument does not 
depend on model Hamiltonian. It should be noted that 
the EMC method, as well as other extended ensemble 
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FIG. 7. Finite-size scaling plot of the domain-wall area in 
4d± J Ising SG model after subtraction of the leading correc- 
tion to the scaling. The critical temperature is used the result 
of the scaling analysis for the domain-wall free energy. The 
exponent v is found to be 0.94(2), consistent with the previous 
estimation. The estimated irrelevant exponent ui = 1.86(77) 
agrees with that obtained from the AFr scaling. The slope 
is estimated 0.94(2), suggesting d s — 3.13(2). 

methods, is also applicable to randomly frustrated spin 
systems, while the cluster-flip based method is restricted 
in non-frustrated models. Another extension would be 
concerned with the choice of the boundary conditions. 
In this paper, we have described the case for the fixed 
spin BC, but it is straightforward to extend it to other 
type of BCs. It is only necessary for boundary conditions 
to be expressed by a countable variable, while the degree 
of freedom of the model system is not restricted. 

We also discuss boundary condition for SG systems. 
Let us comment on related studies. A similar coupled 
replica system has been-studied analytically by mean- 
field variational methodtj, where two replicas are cou- 
pled with each other by fixing the value of overlap be- 
tween surface spins of these replicas. The system stud- 
ied roughly corresponds to the present replica boundary 
model by choosing appropriate parameters. It is pre- 
dicted that an excess free energy due to the effective 
coupling is proportional to L d_5 / 2 , which accidentally co- 
incides with the upper limit of the droplet scaling theory 
in the four dimensional case. Our estimation of the stiff- 
ness exponent is not compatible to that predicted from 
the variational calculation. 

Recently a new boundary condition, called thej-jiaive 
boundary condition, has been proposed in 2D IsingEil and 
XYl3 spin glass models, independently. In these studies, 
they minimize energy of a whole system under the free 
boundary condition. Using the obtained boundary spin 
configuration as a reference system, a twisted boundary 
condition is prepared by flipping the sign of spins on one 
surface. The ground state energy of such system is al- 
ways higher than that of the reference system. They 
claimed that this non-negativity is an evidence of intro- 



ducing correctly a domain wall into the system. It is 
doubtful whether such boundary conditions defined at 
zero temperature is also relevant to the ordering at fi- 
nite temperatures. This is-because many SG systems in- 
cluding both short rangeO and mean field modelso are 
expected to exhibit chaotic nature; namely spin configu- 
rations at finite temperatures differ from those at T = 
in larger scale than the so-called overlap length. Further, 
the replica boundary condition takes an advantage from 
the naive one in a practical sense, because the former 
does not need the ground state calculations. This fact 
makes our investigations easier in three or high dimen- 
sional systems, where the ground states are hardly found 
for suitable large system due to NP hardness. 

The present method has successfully been applied to 
the 4c? ± J Ising SG model under the replica boundary 
conditions. The average of the domain- wall free energy 
AFr over samples, not the variance as used in the stan- 
dard DWRG study, exhibits very clear crossing at the 
critical temperature, implying that it is a good indica- 
tor of the SG transition. It is noted that the replica BC 
is crucial for providing the non negativity of A-Fr. We 
expect that, when the system has a well defined rigidity 
in the ordered phase, the AFr analysis works well even 
in the case where the Binder parameter does not show a 
crossing at T c . In such systems, the short range SG mod- 
els with the field are one of the most attractive problems 
in the SG study. As a byproduct of the RBC, we can 
argue the domain-wall area in the SG phase. We have 
estimated the stiffness exponent 9 and the surface di- 
mension d s of the domain wall in the Ad Ising SG phase 
independently. The latter value lies significantly above 
the trivial surface dimension d — 1, meaning that the do- 
main wall is rough, while both 9 and d s coincide with 
d — 1 in the ferromagnetic Ising models. 

Finally we make a comment on distribution of AFr 
over samples P(AFy{), whose typical results are shown 
in Fig. ^. To our surprise, the distribution functions of 
different sizes, when scaled by their first moment, lie top 
on each others in the SG phase. Another remarkable ob- 
servation is that the scaling function is approximated by 
a Gaussian function; namely it approaches to a nonzero 
value as its argument goes to zero. These results, simi- 
lar to those obsepaefj in 2d and 3d Ising SG models at 
2ratureE3o, are consistent with the droplet 



zero t 



a 



picturet 

The question of whether many equilibrium pure states 
exist or not in the SG phase has still remained contro- 
versial , ..For the system of present interest, some MC 
studiesoC3 have supported the existence of the multi- 
ple pure states, namely the mean field picture, while 
the MigdaJ-Kadanoff approximation for the short range 
SG modeled has claimed that the asymptotic size scale 
to detect the correct thermodynamic properties is far 
from those investigated in the MC simulations. As men- 
tioned in Sec. 



Ill 



the replica BC used in the present 
work prefers a self-overlap configuration in the two repli- 
cas. Correspondingly, under the replica antiperiodic BC, 
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FIG. 8. Scaling plot of the distribution function of the do- 
main-wall free energy with T ~ 1.6. Both axes are scaled by 
the first moment of the distribution. The solid curve is ob- 
tained by fitting the scaling function to a Gaussian formula. 
The raw data is shown in the inset. 

there likely appear such configurations with a domain 
wall which lies in one of the two replicas and separates 
one configuration from its time-reversal one. Therefore, 
our results mentioned above strongly suggest that na- 
ture of low-lying excitations within one pure state is as 
expected in the droplet theory. Our data along, however, 
cannot exclude the possibility that there are many pure 
states. 

In conclusion, we have proposed a MC method which 
enable us to estimate the free-energy difference and have 
successfully applied it to Ad ± J Ising SG model. Our 
value of T c is in good agreement of the previous results 
obtained from the numerical simulations and the series 
expansions. We have presented estimates of two expo- 
nents, the stiffness exponent and the fractal dimension. 
We have also found that low-lying excitations as expected 
in the droplet theory are realized within one pure state in 
the SG phase, though we cannot rule out the possibility 
that there exist many pure states. 
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APPENDIX A: SETTING TEMPERATURE 
POINTS FOR THE EXCHANGE MC METHOD 

In this appendix we propose a practical way to deter- 
mine temperature set which is needed in the exchange 
MC method. For simplicity, we consider a procedure for 
setting a temperature point (3 n between two fixed ones, 
Pn-i and /3„+i. Our criterion is that acceptance proba- 
bilities for the exchange trial with both neighboring tem- 
peratures become equal: 

[fln-l ~ /?„) Wn-l) - E(Pn)) = C, 

(Pn - Wn) - E(P n+1 )) = C, (Al) 

where C and (3 n are unknown constants. A formal solu- 
tion for p n is given by 

Pn = 9(0n) = \~~~pTr T X (&-l^O0»-l) 

E(p n -l) - E(p n+ i) 

-P n+1 E(f3 n+1 ) - S(/3„)(/3„-i - /3„+i)). (A2) 

Regarding 0' = g((3) as a map of (3 to 0', we find an 
fixed point of period 2 with (3 n +i = g(Pn-i) ancl Pn-i — 
g(/3 n +i)- Therefore, we expect a repulsive fixed point 
between and (3 n +i- A new mapping to obtain the 

fixed point is given by 

Pn(t + l) = ^(P n (t)+g((3 n (t))), (A3) 

where t is the iteration step. This iteration scheme can 
be extended straightforwardly to the case for multiple 
temperature points. The whole set of temperature is di- 
vided into two groups with even-n and odd-n. Using the 
iteration scheme, temperature points of the one group 
are updated with the other group fixed, alternatively. In 
actual iterations, the initial temperature points {/3 n } are 
set in a suitable way, for example, equidistant f3. The en- 
ergy E(f3) at the initial set of {3 is roughly estimated by 
short MC simulation and the energy at any temperature 
between [3\ and /3m is assumed to be obtained from the 
MC data, for example, by interpolation technique. The 
convergence of the iteration is rapidly achieved in many 
systems we have investigated. 

From our experiences so far, efficiency of the EMC 
method is rather insensitive for the choice of tempera- 
ture points, when it is applied to systems, such as spin 
glasses, with non-diverging specific heat at the phase 
transition. This fact that it is not necessary to specify 
any parameters before main simulation is, in fact, one of 
big advantages of the EMC method against the other ex- 
tended ensemble methods such as the multicanonical MC 
method and simulated tempering method. Nevertheless 
we emphasize that a little effort on preparing the temper- 
ature points by pre-MC runs following the prescription 
described above ensures the acceptance ratio almost in- 
dependent of temperature and so is quite useful. 
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